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Abstract 

Background: Inflammation triggered by infection or injury is tightly controlled by glucocorticoid hormones which 
signal via a dedicated transcription factor, the Glucocorticoid Receptor (GR), to regulate hundreds of genes. However, 
the hierarchy of transcriptional responses to GR activation and the molecular basis of their oftentimes non-linear 
dynamics are not understood. 

Results: We investigated early glucocorticoid-driven transcriptional events in macrophages, a cell type highly responsive 
to both pro- and anti-inflammatory stimuli. Using whole transcriptome analyses in resting and acutely lipopolysaccharide 
(LPS)-stimulated macrophages, we show that early GR target genes form dense networks with the majority of control 
nodes represented by transcription factors. The expression dynamics of several glucocorticoid-responsive genes are 
consistent with feed forward loops (FFL) and coincide with rapid GR recruitment. Notably, GR binding sites in genes 
encoding members of the KLF transcription factor family colocalize with KLF binding sites. Moreover, our gene 
expression, transcription factor binding and computational data are consistent with the existence of the GR-KLF9-KLF2 
incoherent FFL. Analysis of LPS-down regulated genes revealed striking enrichment in multimerized Zn-fingers- and KRAB 
domain-containing proteins known to bind nucleic acids and repress transcription by propagating heterochromatin. This 
raises an intriguing possibility that an increase in chromatin accessibility in inflammatory macrophages results from broad 
downregulation of negative chromatin remodelers. 

Conclusions: Pro- and anti-inflammatory stimuli alter the expression of a vast array of transcription factors and chromatin 
remodelers. By regulating multiple transcription factors, which propagate the initial hormonal signal, GR acts as a 
coordinating hub in anti-inflammatory responses. As several KLFs promote the anti-inflammatory program in 
macrophages, we propose that GR and KLFs functionally cooperate to curb inflammation. 

Keywords: Transcriptional regulation. Glucocorticoid receptor. Inflammation, Feed forward loops. 
Gene regulatory network, KLF transcription factors 



Background 

Preserving homeostasis is the primary function of the in- 
nate immune system that detects "danger and stranger" 
signals and eliminates invading microorganisms, re- 
sponds to irritation and injury and eventually initiates 
tissue repair. Innate immune cells, such as macrophages, 

* Correspondence: chinenovy@hss.edu; rogatskyi@hss.edu 
^Equal contributors 

^Hospital for Special Surgery, The David Rosensweig Genomics Center, 535 
East 70th Street, New York, NY 10021, USA 

^Graduate Program in Biochemistry, Cell and Molecular Biology, Weill Cornell 
Graduate School of Medical Sciences, 1300 York Avenue, New York, NY 
10021, USA 

Full list of author information is available at the end of the article 

(3 BioMed Central 



neutrophils and dendritic cells constantly sample their 
environment for lipopolysaccharides (LPS), double and 
single-stranded nucleic acids, microbial proteins and 
other broad molecular patterns that are not normally 
present in eukaryotes and, in response, produce cyto- 
kines and chemokines that attract additional immune 
cells to the site of infection or injury. Normally a pro- 
tective response, excessive or persistent inflammation is 
associated with tissue damage and needs to be regulated. 
Indeed, numerous mechanisms have evolved that control 
inflammation at multiple levels. Systemically, inflam- 
matory stimuli activate neuro-endocrine circuitry that 
triggers the production of glucocorticoids by the adrenal 
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glands, which ultimately attenuate the expression of in- 
flammatory cytokines [1]. These potent anti- inflammatory 
properties made glucocorticoids exceedingly common for 
managing a wide range of autoimmune and inflammatory 
conditions such as rheumatoid arthritis, systemic lupus 
erythematosus, inflammatory bowel disease, psoriasis and 
multiple sclerosis [2]. 

Glucocorticoids signal through the glucocorticoid recep- 
tor (GR) - a ubiquitously expressed ligand- dependent 
transcription factor (TF) of the nuclear receptor (NR) 
superfamily. GR regulates transcription by either binding 
directly to specific DNA sequences known as glucocortic- 
oid response elements (GREs) or by tethering to other 
DNA-bound regulators, such as Activator Protein (AP)1, 
Nuclear Factor (NF) kB and Signal Transducers and 
Activators of Transcription (STAT) family members [3]. 
Although GR is expressed in all immune cells, the phy- 
siological outcomes of GR activation are highly ceU type- 
specific: for example, glucocorticoids are anti-apoptotic in 
neutrophils, but pro-apoptotic in eosinophils, dendritic 
and some T-cells [4]. Moreover, prolonged glucocorticoid 
treatment induces cells polarization toward a new pheno- 
type with either pro- or anti-inflammatory properties [5,6]. 

The analyses of glucocorticoid-regulated transcriptomes 
paint a complicated picture encompassing hundreds of 
up- and down-regulated genes that vary in different cell 
types and populations and over time of glucocorticoid ex- 
posure [7,8]. Although shared GR target genes certainly 
exist, system-specific regulators and pathways drastically 
affect transcriptional outcomes, response dynamics and 
relative activities of such shared genes and their products. 
The existence of intricate inter-protein and inter-pathway 
interactions contributing to the NR-mediated gene regula- 
tion has been proposed almost 30 years ago [9]. The struc- 
tural analysis of NR transcriptional networks, however, 
was lagging due to the lack of genome-wide data and li- 
mited availability of analytical tools. More recently, studies 
in bacteria and yeast have defined specific patterns of 
functional interactions ("network motifs") between inter- 
dependent TFs and provided a computational framework 
for the analysis of gene expression data to identify such 
motifs [10-12]. The simplest motif - a positive or negative 
auto- regulatory loop - consists of a single TF that regu- 
lates its own expression. More complicated feed forward 
loops (FFL) involve three factors: a signal- responsive mas- 
ter regulator TF, an intermediate TF controlled by the 
master TF and a jointly regulated gene under the control 
of both the master and intermediate TFs [12]. Depending 
on the specific activities of the master and intermediate 
TFs (activation or repression) and the response thresholds 
of participating genes, the dynamics of the FFL transcrip- 
tional outputs vary, yielding unique expression patterns 
for various TF and target gene combinations. Dynamic 
responses elicited by FFLs deviate from simple gene 



regulation providing for unusual control mechanisms 
that are responsible for noise filtering, fold-change 
sensing, pulse generation and transcription response 
acceleration [11,13,14]. 

Gene expression networks can be represented as graphs 
with TFs and other expression regulators acting as nodes 
and functional interactions between regulators and bet- 
ween regulators and their targets as directional edges 
[12,15]. Of particular interest are the networks and net- 
work motifs in which TFs regulate each other, or them- 
selves act as transcriptome organizers and ultimately 
determine the topology of the entire network [12]. Several 
computational algorithms and experimental approaches 
have been successfully applied to map global trans- 
criptional networks and identif)^ novel functional motifs in 
various organisms [16,17]. In metazoans, this analysis is 
often complicated by the lack of information about the 
edge identity (not all targets for a given TF are known, 
some known "targets" are not regulated directly) and di- 
rection (a TF can either activate or repress the same gene 
in a tissue-specific manner). To complicate matters 
further, the role of intermediate TFs can be fulfilled by miR- 
NAs or regulators of RNA translation and stability [18,19]. 
Thus, dissecting regulatory networks requires a combi- 
nation of computational and experimental approaches. 

We reasoned that a highly branched response to gluco- 
corticoids is determined by the early transcriptional 
events. Here, we focused on the regulatory network elic- 
ited by an acute stimulation of mouse macrophages with 
glucocorticoids and/or LPS. Combined with high-reso- 
lution kinetic experiments and dynamics modeling, this 
analysis enabled us to dissect early post-stimulation events 
prior to extensive signal propagation, which usually masks 
the bona fide response to GR activation by a web of 
secondary effects. 

Results 

Transcriptome analysis of mouse macrophages exposed 
to acute glucocorticoid and LPS stimulation 

To analyze early regulatory events initiated by glucocorti- 
coids and inflammatory stimuli we treated BMMO with 
either ethanol vehicle (U), LPS (L), Dex (D), or a combin- 
ation of the two (L + D) for 1 h, isolated and sequenced 
PolyA-enriched RNA as described in (Additional file 1). 
The sequencing results are summarized in Additional 
file 2: Table SI. 

To uncover the regulatory patterns in gene expression 
data, we performed /c-mean cluster analysis of ANOVA- 
filtered differentially expressed genes (see Additional file 1). 
To minimize magnitude-based clustering, the log2-trans- 
formed expression values were first converted to Z-scores 
as in Z = (Xgem-X) / a where Xggne is an expression value 
for a given gene at a given condition, X is an average 
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expression value across conditions and a is a standard 
deviation of X across conditions. The optimal number of 
clusters was determined using the "elbow" method by 
plotting within-cluster variance vs. the number of clusters. 
The analysis was performed using Euclidian distance with 
progressively increasing number of clusters from 8 to 14 to 
determine a stable configuration using the cluster analysis 
module of STATISTICA 8.0. We grouped ANOVA-filtered 
data into 12 clusters of co-regulated genes (Additional file 
2: Table S2). We further evaluated the significance of differ- 
ences in gene expression within clusters by the Mann- 
Whitney test and validated the results for a limited number 
of genes by RT-qPCR. As many glucocorticoid- and LPS- 
responsive genes have been previously characterized by us 
and others [20-22], we selected for validation either poorly 
characterized or novel target genes. Based on patterns of 
co-regulation, we grouped these clusters into four larger 
categories (Figure 1). 

I) Genes activated by either LPS (Cluster 2) or 
glucocorticoids (Cluster 7) 

The majority of genes in these two clusters are upregulated 
by either LPS or Dex independently, with little evidence 
for inter- treatment interactions at 1 h. The LPS -induced 
Cluster 2 (Figure 1, Mann- Whitney, Pu-l = 4.11*10"^^) 
contains genes encoding pro- and anti-inflammatory cyto- 
kines and chemokines (IllO, CxcU, 3, 5 and 7, Ccl7 and 
Tnfsf9), TFs involved in stress response (Maff, Ets2, Fosl2 
and Kdm6b) and proteins involved in TLR signaling (Tlr2, 
Cdl4 and Cd40) and signal transduction (Itpkc, Rabgefl, 
GbpSa). 

Cluster 7 contains Dex-induced genes (Pu-d = 0.0013), 
including several well-characterized GR targets such as 
TFs Perl and Klf9, immunophilin FkbpS, potassium chan- 
nel Kcnk6. In addition, this cluster includes several genes 
whose regulation by Dex has not been previously reported: 
Interleukin 15 receptor alpha (lUSra), the Wnt pathway 
receptor Fzd4, the TF Klf2 and chemokine Cell 7. 

II) Genes co-activated by LPS and Dex 

These genes display either predominantly additive {cluster 
1) or synergistic {cluster 6) activation by LPS and Dex. 
Several genes in these clusters are previously characterized 
GR targets including Duspl, Nfil3 and Cited2 {cluster 1) 
and Mt2 and Pfkfb3 {cluster 6), whereas the glucocor- 
ticoid and LPS responsiveness of several others, such as 
histone chaperon Jdp2 {cluster i), has not been reported 
previously. 

III) Genes induced by LPS and repressed by Dex 

This large group of genes is represented by clusters 3, 4 and 
5. Cluster 3 contains LPS-induced genes (Pu-L = i-95no-^') 
expressed at relatively high level in resting BMMO. The 



basal expression of these genes is significantly more sen- 
sitive to hormonal treatment (Pu-d = 0.0107) than their 
LPS-induced expression. This cluster encompasses a 
number of inflammatory cytokines (Ccl2, 3 and 4, Tnf, 
Tnfaip2), TFs (lerS, Junb, Bcl6, Prdml and Irfl) and pro- 
teins involved in signal transduction (Gadd45b, DuspS, 
Rasgeflb). Interestingly, several genes in this cluster (Ccl2, 
3 and 4, Tnf) are characterized by the presence of the 
stalled RNA Pol II near the transcription start site in unin- 
duced conditions and are activated primarily at the level 
of the Pol II pause release during early elongation [23-25] . 
Cluster 4 combines a heterogeneous group of genes with 
low basal expression (Puduster 4<uciuster3 = 0.0003) that are 
strongly induced by LPS and includes inflammatory cyto- 
kines (nia, nib, CxcllO, Ifnbl, Tnfsf4 and Illf9) and other 
direct mediators of inflammation (Ptgs2, mIR-155 host 
gene, Hspla), TFs (Egr3) and signaling molecules (Gpr84, 
Areg). Cluster 5 contains many genes whose LPS induction 
is strongly inhibited by Dex treatment (Pl-(l+d) = 0,02) in- 
cluding several cytokines and chemokines (1112b, Lif, lUrn 
and I117ra ) and TFs (Mxdl, Etv3, Klf7 and Etsl). 

Cluster 8 contains statistically heterogeneous previously 
reported (Atf3, Egr2 and Ier3) as well as novel (End and 
Bhlhe40) targets for GR-mediated repression that are 
largely unaffected by LPS treatment. Thus, this cluster is 
formally outside of group III. 

IV) Genes downregulated by LPS 

LPS -repressed genes were separated into 3 clusters based 
on the combined effect of Dex and LPS on gene ex- 
pression. The LPS-mediated downregulation is either 
weakly potentiated {cluster 12) or antagonized {cluster 11) 
by Dex. The functions of the majority of these genes are 
poorly understood or unknown, however, 32% of genes in 
cluster 11 and 20% in cluster 12 encode uncharacterized 
C2H2 Zn-finger proteins implicated in transcription and 
chromosome maintenance. 

Expression of genes in cluster 10 is downregulated by 
Dex and LPS in an additive manner. At least one gene in 
this cluster (Angptl4) is a known GR target. Several genes 
encode regulators of immune cells activities including Ritl 
and Cd300lb. Cluster 9 contains Dex-induced genes that 
are weakly repressed by LPS including previously reported 
Ddit4, Arl4d and Sikl. 

Cluster validation by RT-qPCR 

Two representative genes in each cluster were chosen 
for validation. Total RNA was isolated from treated (D, 
L, or L + D for 1 h) and control BMMO and transcript 
levels for indicated genes (Figure 1) were determined 
by RT-qPCR using Actl or Hprt as housekeeping control 
genes. As several LPS-induced, Dex-repressed genes 
coding for various cytokines that we found in clusters 3 
(Tnf, Ccl2, 3 and 4), 4 (flla and b) and 5 (Lif, Niacrl, 



Chinenov et al. BMC Genomics 2014, 15:656 
http://www.bionnedcentral.conn/1471 -21 64/1 5/656 



Page 4 of 1 9 



I Cluster 2, N=89 

■ (Tnfaip3, CxcH, 3) 



Cluster 7, N=32 

(Perl, Fkbp5, Tsc22d3) 



II Cluster 1, N=22 

" (Nfil3,Jdp2) 



Cluster 6, N=23 

(Pfkbf3, S100a10) 





Cluster 3, N=43 

(Tnf, Ccl2, Ccl3, Ccl4) 



Cluster 4, N=89 

(111 b, CxcHO, Egr3) 



Clusters, N=43 

(Lif, Niacrl, Hspalb) 



Cluster 8, N=20 

(Egr2, Bcl3) 




UDLD UDLD 
L 

Etv3 



1112b 

1 



UDLD UDLD 



Atf3 Gadd45a 



1 




\y Cluster 10, N=35 



Cluster 11, N=51 



Cluster 12, N=69 



DID DID 



Cluster 9, N=27 



Angptl4 



UDLD UDLD UDLD ^^UDLD UDLD 

L L L L L 

Cbx8 Zfp623 

m 





Figure 1 Dex- and LPS-regulated genes in BMMO form distinct clusters depending on specific patterns of expression. ANOVA-filtered 
RNA-seq expression values (RPKM) were Z score-transformed and subjected to k-mean clustering. For each cluster, the upper left panel shows the 
cluster average Z score-transformed log2 RPKM for each treatment condition where the central square represents standardized cluster mean, the 
rectangle is the mean +/- standard deviation (SD) and the whiskers are a 95% confidence interval; the upper right panel shows log-transformed raw 
cluster mean expression values. Symbols representing treatment conditions are black - untreated (U), blue - Dex-treated (D), pink - LPS-treated (L) and 
purple - co-treated (DL) cells. The expression of two representative genes per cluster were determined by RT-qPCR with gene-specific primers and 
shown at the bottom. The statistical significance of differences (Mann-Whitney test) is indicated by asterisks as following: - p < 0.05, '''' - p < 0.01, 
p < 0.001 and ns is non-significant. 



Mmpl3) have been previously characterized by us and RT-qPCR closely resembled those determined by RNA- 
others [20-22], we focused on uncharacterized LPS/ seq, with the exception of several weakly expressed 
Dex targets. The expression patterns determined by genes that were not detectably repressed by Dex 
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following a 1-h treatment {e.g., Bcl2, Med21; data 
not shown). 

Glucocorticoid-regulated genes form a highly 
interconnected association network with distinct 
response-specific modules 

To determine the prevalent functional patterns in groups 
of co-regulated genes, we performed gene enrichment 
analysis and visualization using GeneMANIA plugin [26] 
for Cytoscape 2.8 and Exploratory Gene Association 
Networks (EG AN) software [27]. 

Using a list of genes that were up- {clusters 1,6, 7 and 9) 
and down- {clusters 3, 4, 5, 8 and 10) regulated by Dex, we 
generated a consensus association network consisting of 
333 nodes and 8296 edges. To discern underlying data 
structure in the Dex-regulated network, we decomposed 
this network using the Newman-Girvan community clus- 
tering algorithm [28], a divisive procedure that iteratively 
removes network edges with largest "edge betweenness", 
recalculates this metric for a novel network and repeats 
the procedure until the network is split into several 
groups. The Newman-Girvan algorithm disregards edge 
weights and uses only connectivity to define communities. 
Community clustering partitioned Dex-regulated network 
into three unequal modules that were significantly en- 
riched with Dex-repressed genes (Module 1, Figure 2A), 
Dex-induced genes (Module 2, = 18.33, p = 0.0001, 
Figure 2B) or LPS -repressed genes (Module 3, = 15.347, 
p = 0.00046, Figure 2C). 

Network topology analysis of Dex-responsive modules 
revealed significantly greater network densities and clus- 
tering coefficients in Modules 1 and 2 compared to net- 
works generated from the same number of non-expressed 
genes extracted from the same BMMO experiments 
(Figure 2D). Similarly, the average neighborhood connec- 
tivities and average clustering coefficients for Modules 1 
and 2 were considerably greater than the values for non- 
expressors (Figure 2E), indicating that nodes in these mo- 
dules form tight interconnected local groups. A broader 
shared neighbors distribution in Modules 1 and 2 indicates 
high prevalence of shared nodes, suggesting an enrich- 
ment for multiple input motifs. Interestingly, Module 2 
(enriched with Dex-induced genes) was more structured 
than Module 1 (predominantly Dex-repressed genes) 
as evidenced by consistently higher average neighbor- 
hood connectivities and average clustering coefficients 
(Figure 2E). Although Module 3 has a non-random com- 
position, with the exception of network heterogeneity, all 
other analyzed topological metrics were typical of net- 
works composed of randomly selected non-expressed 
genes (Figure 2D). Therefore, we focused the rest of the 
analysis of Modules 1 and 2. Relatively high heterogeneity 
for both modules (0.506 and 0.581; Figure 2D) indicates 
the presence of network hubs - nodes with high degree of 



connectivity. Indeed, 10 most connected nodes in Mo- 
dules 1 and 2 account for 36 and 32% of all edges, res- 
pectively (Figure 2F). Interestingly, eight out of 10 most 
connected nodes in Module 2 are sequence-specific DNA- 
binding TFs (bold in Figure 2F). 

Glucocorticoid response-specific modules are functionally 
distinct 

Using gene ontology (GO) analysis to determine enriched 
gene categories in subsets of functionally related genes, we 
identified 425 enriched GO categories for Module 1, 285 
for Module 2 (FDR corrected p < 10"^) and 77 for Module 
3 (FDR corrected p < 10"^). Only 115 GO categories over- 
lapped between Modules 1 and 2 (Figure 3A, 3C). To fa- 
cilitate visualization and interpretation of these results and 
compare enriched functional categories among groups of 
Dex-regulated genes, we generated GO terms similarity 
networks using Gene Set Enrichment Mapping Cytoscape 
plug-in. Multiple GO categories related to regulation of 
metabolic processes, embryonic and post-embryonic de- 
velopment and regulation of apoptosis and signaling are 
enriched in Module 2 (Figure 3A) that contains a large 
number of Dex-upregulated genes. Notably, 32/285 GO 
categories enriched in Module 2 were related to regulation 
of gene expression, regulation of transcription, sequence- 
specific DNA binding transcription factors. For example, 
negative regulation of gene expression (GO: 10629), ne- 
gative regulation of transcription (GO: 16481), negative 
regulation of transcription - DNA-dependent (GO:45892), 
sequence-specific DNA binding (GO:43564), negative re- 
gulation of transcription from RNA polymerase II pro- 
moter (GO: 122) and transcription regulator activity (GO: 
30528) were all enriched in Module 2, but not Module 1 
(Figure 3C). Overrepresentation of genes coding for re- 
gulators of gene expression in the early Dex-responsive 
transcriptome suggests that GR initiates a transcriptional 
program that relies on the step-wise activation of multiple 
TFs. Only a few categories related to immune/inflamma- 
tory responses have been found in Module 2 (Figure 3A). 

Conversely, the majority of enriched GO categories in 
Module 1, which contains predominantly Dex-repressed 
genes, are related to immune and inflammatory responses, 
signaling and regulation of signal transduction and meta- 
bolic regulation including immune response (GO:6955), 
immune system process (GO:2376), inflammatory re- 
sponse (GO:6954) and regulation of cytokine production 
(GO: 18 17) (Figure 3B, 3C). The fraction of gene expres- 
sion-related GO categories in this module is significantly 
smaller (22/425, x^ = 8.05, df = 1, p = 0.00455) than in 
Module 2. 

Dex-responsive transcription regulators 

We identified 37 Dex-responsive genes whose products 
are involved in the regulation of gene expression (Figure 4). 
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Figure 2 (See legend on next page.) 
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(See figure on previous page.) 

Figure 2 Hormone-regulated genes in BMMO participate in highly interconnected networks with TF control nodes. (A) Ranl<ed list of 
Dex-responsive genes was used as an input for GeneMania networl< building algorithm. The resulting combined network was partitioned using 
the Newman-Girvan algorithm into three Modules enriched with (A) Dex-repressed, (B) Dex-induced and (C) LPS-repressed genes. To simplify the 
view, the edges representing co-expression, co-localization and shared protein domains were hidden. The size of the node is proportional to the 
node connectivity. The nodes are colored according to the magnitude of Dex effect on respective gene expression using MulticoloredNode Cytoscape 
plugin. The nodes predicted by GeneMania are colored gray. Both (D) network-wide and (E) node-specific network topological parameters in Modules 
1 and 2 demonstrate high level of interconnectivity with a large number of shared nodes. (F) Nodes with the highest degree of connectivity in the 
Modules 1 and 2. The genes coding for TPs are shown in bold. 



12 of these genes including TFs Klf2, 4 and 9, Perl, Jdp2, 
Cited2, Nfil3 and Tiparp are upregulated by Dex; 12 
others including TFs Junb, Atf3, Tgifl, Irfl and Bcl3 are 
downregulated; for 13 regulatory proteins {e.g,, Zfpl31, 
Zfp36, Nr4a3, Rela, Nfl<b2, Klf7 and Etsl), the inhibitory 
effect of Dex is apparent only in LPS -induced MO. For all 
Dex-induced and a subset of Dex-repressed TFs, we have 
independently confirmed RNA-seq data by RT-qPCR 
(Additional file 3: Figure SI; also see Figure 1 for Cited2, 
Irfl, Etv3 and Atf3). 

To uncover potential functional interactions between 
Dex-regulated TFs, we treated BMMO with Dex for up 
to 9 h and determined expression levels of a subset of 
Dex-regulated genes by RT-qPCR. We observed a stri- 
king difference in expression patterns over time. Nfil3, 
Cited2, Jdp2 and Perl (Figure 5A) are characterized by 
an accelerated burst phase, with the mRNA level rea- 
ching maximum within 30-60 min and then remaining 
constant (Nfil3 and Cited2) or slowly declining over time 
(Perl and NcoaS). Klf4, Klf9, Tsc22d3 and Ddit4 are 
strongly induced within the first 2 h, and their RNA 
levels continue to increase for the next 6 h (Figure 5A). 
The expression profile of FkbpS exhibits an initial delay, 
followed by a robust and sustained induction (Figure 5A). 
Conversely, Klf2 and Tiparp displayed pulse-like rapid 
upregulation within 1-3 h followed by a decline in tran- 
script level, which in the case of Klf2 reaches baseline 
(Figure 5A); a similar biphasic pattern of expression was 
observed for Tgfb3, I115ra and Mt2 (Figure 5A). Interes- 
tingly, Bcl3, Junb and Tgifl responded with rapid pulse- 
like downregulation followed by a slow return to basal 
expression level, whereas Atf3 was rapidly downre- 
gulated within the first hour and remained repressed 
throughout the time course (Figure 5A). Unexpectedly, 
the Pparg expression was only modestly induced by Dex 
at the early time points, then decreased dramatically by 
3 h and remained low for up to 9 h (Figure 5A). 

The dynamics of expression for several Dex-regulated 
TFs suggests that they are under combinatorial controls 
that involve GR and additional GR targets which either 
cooperate with or antagonize GR actions. As such a model 
implies transcription/protein production of these putative 
GR targets, we first examined the expression of Dex- 
regulated genes in the presence of a protein synthesis 



inhibitor cycloheximide (Chx). Treatment with Chx up for 
to 3 h had no dramatic effects on Dex-mediated re- 
gulation of Klf9, Tcs22d3, Tgfb3 and Bcl3 (Figure 5B), 
suggesting a direct regulation by GR that does not rely on 
synthesis of additional proteins. Conversely, the expres- 
sion of Atf3 became refractory to Dex in the presence of 
Chx (Figure 5B) suggesting that additional proteins in- 
duced by Dex rather than GR itself are likely to directly 
regulate this gene. For several Dex-responsive genes (Klf2, 
Klf4, Nfil3, Tiparp, Tgifl), however, treatment with Chx 
dramatically upregulated their basal expression, complicat- 
ing the assessment of the relative contribution of direct vs. 
indirect effects of GR to their regulation (Additional file 3: 
Figure S2) and necessitating an alternative approach. 

Temporal dynamics of hormone-regulated gene 
expression is consistent with feed forward logic 

The dynamics of the transcriptional response of several 
genes to Dex imply the existence of some feedback me- 
chanism that limits activation by GR yet, at the same time, 
is GR-dependent. Because Chx elicits many off-target ef- 
fects and does not enable us to discriminate between the 
secondary targets of GR and those jointly controlled by 
GR and a GR-regulated TF, we performed dynamic mo- 
deling of expression data in an attempt to discern specific 
regulatory patterns. Several mechanisms including positive 
and negative autoregulation, positive and negative feed- 
back and feed forward loops (FFL) could account for de- 
viations from a simple model with a single TF regulating 
gene expression via a single DNA binding site [12]. The 
kinetics of GR expression following Dex stimulation is not 
consistent with auto-regulatory models (Figure 5A). De- 
pending on the overall regulatory outputs and activities of 
individual FFL components, two types of FFL have been 
recognized - coherent (C-FFL) and incoherent (I-FFL). 
In type 1 I-FFL (Il-FFL), the activating master TF and a 
repressing intermediate regulator have opposite effect on 
a jointly regulated gene. Dynamic modeling and experi- 
mental studies of Il-FFL dynamics have demonstrated 
several properties of this network motif, including its 
ability to produce sharp pulse-like activation of a jointly 
regulated gene with a fast relaxation time [11,13]. Several 
GR-regulated genes, including Klf2, Tiparp, Tgfb3, Mt2 
exhibit pulse-like kinetics at constant Dex exposure 
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(See figure on previous page.) 

Figure 3 Distinct functional GO categories are enriched in Dex-regulated networic Modules in BMMO. GO categories enriclied among 
Dex-regulated genes in Modules 2 (A) and 1 (B) (p < 10"^ liypergeometric test) were used to create GO terms similarity networks where the 
nodes represent GO gene sets and the edge length between two nodes is proportional to the fraction of shared genes in categories. Individual 
clusters were manually repositioned to simplify network layout and colored as indicated; the broad top GO terms are colored grey. The word tags 
were generated with WordTag Cytoscape plugin [29] to represent frequency-based semantic summary of words in GO category titles found in 
individual clusters and the size of the nodes is proportional to the number of genes in the respective GO gene set. (C) GO enrichment analysis of 
individual GO categories in Modules 1 and 2 relative to the mouse genome. The FDR-corrected p-values (hypergeometric test) are shown for 
each bar. 



(Figure 5A). Near-baseline relaxation of Klf2 expression 
suggests that this gene is under control of both GR and a 
strong Dex-induced repressor. Because GR is largely in- 
active in the absence of ligand, glucocorticoids act as a low- 
latency on-off switch eliminating the need to correct for a 
baseline activity of GR. The dynamics of Klf2 and repressor 
(R) expression is described by the ordinary differential 
equations (2) and (3) (see Additional file 1) [11,13]. 

Assuming equal degradation rates of Klf2 and R, these 
equations can be solved analytically (equation (1) and 
[13]) and used to fit the expression data for Klf2. When 
limited to the early data points (up to 4 h), the expres- 
sion data fit very well to the predicted expression pattern 
(Figure 5C, = 0.9225), however, at the later time, when 
the contribution of degradation rates becomes sig- 
nificant, the quality of fit decreases (R^ = 0.5491). Using 
parameter estimates derived from equation (1) fitting, 
we solved equations (2) and (3) numerically. Figure 5C 
shows a good concordance between the experimental 
data for Klf2 expression as determined by RT-qPCR and 
the numerical solution {R^ = 0.8317) that describes the 
dynamics of the Il-FFL. This result strongly suggests 
that Klf2 and, potentially, several other GR targets that 
exhibit similar expression dynamics (Tiparp, Tgfb3 and 
Mt2) are jointly regulated by GR and a GR-induced re- 
pressor via the Il-FFL network motifs. 

Numerical solutions of equations 2 and 3 also provide 
a theoretical prediction of an intermediate repressor dy- 
namics in the GR-R-Klf2 I-FFL. Interestingly, among 
several known transcription repressors activated by GR, 
the expression kinetics of Klf9 fits the best to the pre- 
dicted model (Additional file 3: Figure S3). The dynam- 
ics of the GR-R-KLF2 FFL can be tested by perturbing 
the concentration of a hypothetical intermediate repres- 
sor which should uncouple the FFL thus shifting peak- 
like FFL-mediated kinetics to simple monotonous kinetics 
eventually converging to a steady-state level. To test the 
role of Klf9 as a potential GR-activated repressor of 
Klf2 transcription, we derived M from Klf9-KO mice 
[30] kindly provided by Dr. Simmen and treated them 
with 100 nM Dex as above. Interestingly, in Klf9 null 
BMMO the peak-like Klf2 induction profile "degenerated" 
to monotonous activation kinetics that plateaued at a 
steady state level by 3 h (Figure 5D), replicating the 



profiles of previously reported uncoupled experimental I- 
FFLs [11,31], consistent with the proposed role of KLF9 as 
an intermediate repressor in the GR-KLF9-KLF2 I-FFL. At 
the same time, deletion of Klf9 did not affect the expres- 
sion dynamics of either GR itself or the GR target gene 
Tsc22d3 with a simple monotonous activation profile 
(Figure 5D). 

GR is recruited to binding sites associated with 
Dex-regulated genes 

The FFL gene regulatory circuitry predicts that the master 
TF binds DNA to regulate transcription of both FFL 
nodes. Using several published mouse ChlP-seq datasets 
of acute GR recruitment [32,33] we interrogated Dex- 
induced genes for the presence of GR binding sites within 
the gene and 15 Kb upstream and downstream of the gene. 
With the exception of Sox4 and Klf4, all genes encoding 
Dex-induced TFs contained at least one peak within the 
analyzed intervals (Figure 6A, Figure 4). To compare these 
peaks to those in MO, we analyzed GR recruitment by 
ChlP-seq using untreated MO as a control and identified 
16,657 peaks induced by a 40-min Dex exposure at 2% 
FDR. Selective comparison of binding site distributions re- 
vealed a high level of concordance between Dex-induced 
peaks in MO and those previously described in adipocytes 
[32] (Figures 6A, 7B and Additional file 3: Figure S4) and a 
partial overlap with a GR cistrome in MO polarized with 
high dose long-term glucocorticoid exposure [34]. By 
ChlP-qPCR, we detected GR recruitment as early as 
40 min post Dex treatment at multiple putative GR bin- 
ding sites, including those at Perl, Cited2, Klf2, Klf9, Nfil3, 
Jdp2, Tiparp and Ncoa5 (Figure 6B). These observations 
correlate well with the expression data (Figures 1 and 5A). 
Although Klf4 was strongly induced by Dex, no gluco- 
corticoid response elements (GREs) near the gene has 
been previously reported. We performed a scanning ChIP 
with evenly spaced primers within the Klf4 gene and 
several primers flanking potential GR binding sites 
(Figure 6C). Two of the putative GREs located at -3830 bp 
(gGcACAgcaTGTaTC) and +5896 bp (aGaACAgaaTG 
Tagttc) relative to the Klf4 transcription start site recruited 
GR following a 40-min treatment with Dex (Figure 6C), 
consistent with the notion that, similar to genes shown in 
Figure 6B, GR is likely to regulate Klf4 directly. 



Chinenov et al. BMC Genomics 2014, 15:656 
http://www.bionnedcentral.conn/1471 -21 64/1 5/656 



Page 10 of 19 



Feature 


U 


D 


L 


(L+D) 


Associated 


ID 




lu^^y rxr rvivi / 




luy^y r\r rvivi j 


GRE 


l\IT4 


1 /I Q 


4.32 


i.by 


3 "7Q 

o. /y 


V 

Y 


reri 


Z.z>4 


■ ^-^^ 


T Q"7 

z.y / 


C 2/1 

D,o4 


Y 

T 


Litedz 


D.D/ 


1^ o.ii 


-7 r 
/.D 


Q C/l 

y.D4 


V 
Y 


Kity 


1.51 


3.3 


1.28 


2.87 


Y 


l\ITz 


o./b 


■ D.4y 


/I C 1 


C AC 

D.4b 


Y 


M-FilQ 

NTIIo 


T "71 
Z./l 


HI 4. ZD 


4.4/ 


D.yb 


v 
Y 


Tiparp 


Z.DZ 


4.Uj 


4.0 J 


c ns 
D.Uo 


V 
Y 


Jupz 


5.DD 


4.bD 


b.zy 


/.Zo 


V 
Y 


Sox4 


U.D J 


1 f^7 

1.67 


U.4i7 


n f^A 

U.D4 




NcoaS 


z.yo 


i.o/ 


i.lz 


Q OQ 
D.OD 


v 
Y 




o.yi 


■v 4.D/ 


D.ZD 


C QT 
D.bZ 


V 

Y 


Mr/1 o1 

IN r 4ai 


Z.o/ 


d.Zd 


C 73 


C 37 
b.DZ 


V 

Y 


Kei 


Z.O 


Z.DZ 


3.Z 


/I Q3 


V 

Y 


I\ITD 




D.ZZ 


"7 CQ 
/.DO 


"7 3 
/.D 


V 

Y 


Z.Tpiol 


l.D 




3 71 


O 
J 




Z.Tpob 


b.l 


C 1 1 

b.ll 


Q Q 3 
y.DO 


Q no 
y.Uo 


V 

Y 


Nr4a3 


0.17 


0.13 


1.07 


0.4 




Rela 


D.O/ 


C 17 
J.Z / 


C /I 3 

b.4o 


C 1 1 

b.ll 




NTKDZ 


3.14 


4,yD 


C 71 

b. /z 


C 3Q 




Myc 


Z.Z 


z.Ui 


d.Ud 


A C 

4.b 




l/lf"7 

l\IT/ 


Z 


1 o 
l.o 


3 C"7 

d.d/ 


T "7C 

z./b 




tisi 


1.34 


i.zy 


Z.DO 


1 Q7 
l.b/ 




bati 


4.ZO 


3 Q"7 

o.y / 


A QC 

4.yb 


/I QC 




Qc1 SI1 


l.cS4 


i.Db 


1 


n Q3 




USTi 


3.U1 


/I "7 

4. / 


/111 

4.11 


3 QC 




ivixai 


z.yb 


z.oy 


C /I Q 

D.4o 


/111 
4.11 




bgrl 


z.o4 


1 QQ 
l.OO 


b.y 


C CI 

b.bl 


V 

Y 


Dnine4u 


T 1 Q 


l.Z) 


z.yz 


7 1 C 

z.lb 


V 

Y 


Bcl6 


^ 1Q 


IB 9 ZLQ 


4.58 


4,17 


Y 


Etv3 


3.49 


2.73 


5.08 


4.13 




Sertadl 


4.35 


3.37 


4.81 


4.77 




Tgifl 


4.19 


3.14 


5.56 


5.12 




SertadS 


3.23 


^ 2.1 


1.72 


1.3 




Irfl 


4.66 


■ 3.52 


7.48 


6.46 




Bcl3 


4.6 


I 3.35 


6.07 


5.5 




AtfB 


6.78 


^ 5.48 


8.06 


7.5 




Junb 


5.88 


^^.51 


9.23 


8.67 





Dex-induced genes are shown in red, Dex-inhibited - in green. For genes shaded grey, the 
inhibitory effect of Dex is apparent only following LPS induction. 

Figure 4 Dex-responsive regulators of gene expression. 



We then used a set of Dex/LPS-regulated genes to 
analyze the distribution and enrichment of binding sites 
for TFs that were present in Chip-Enrichment Analysis 
(ChEA) database [35]. The ChEA database contains cu- 
rated published genome-wide datasets of TP binding 
sites in human, mouse and rat. After filtering out TFs 
that were not expressed in MO (RPKM < 1), we noted 
that binding sites for several Dex-responsive TFs, such 
as KLF2, KLF4, ATF3, EGRl, CEBPp and IRFl are 
enriched among Dex/LPS-regulated genes. Interestingly, 
binding sites for PPARy, whose expression was inhibited 



upon prolonged Dex treatment, were found near the 
majority (22/37) of Dex-responsive gene expression 
regulators (Figure 4) and highly enriched among Dex/ 
LPS -regulated genes in general. 

We then determined the frequency of genes associated 
with binding sites for several TFs identified by ChEA and 
induced by Dex in individual clusters of Dex/LPS-regu- 
lated genes (Figure 1). We defined a binding site as 
'gene-associated' if its genomic intervals overlapped with 
genomic intervals encompassing all mouse genes anno- 
tated in mm9 +/- 15 Kb by at least one nucleotide. In 



Chinenov et al. BMC Genomics 2014, 15:656 
http://www.bionnedcentral.conn/1471 -21 64/1 5/656 



Page 11 of 1 9 



0123456789 



0123456789 



u 



0123456789 
c Tsc22d3 



0123456789 
Perl 



0123456789 



Ddit4 



dita 



0123456789 



Jdp2 



0123456789 



FkbpS 



0123456789 



^ ♦ ♦ 



0123456789 



Tiparp 



0123456789 



0123456789 



0123456789 



Tgfb3 



0123456789 



0123456789 



0123456789 



S 2.0 

Q. 1.5 



Tgifl 



0123456789 



g 0.6 
aj 0.4 
0.2 



0123456789 



Pparg 



0123456789 



0123456789 

Dex, hours 



0123456789 



0123456789 

Dex, hours 




□ Control 
Dex Ih 
m Dex 2h 
■ Dex 3h 



GR- 

L 



|Klf2 



• Klf2 

fitted solution (/?' = 0.5491) 

— fitted solution (/?' = 0.9225) 

— numerical solution {R^ = 0.8317) 




Dex, hours 




Nr3cl 



3 4 5 6 

Dex, hours 



2 3 4 5 6 

Dex, hours 



Figure 5 The dynamics of hormone-responsive gene expression is consistent with FFL networic motifs. (A) BMMO were treated with 1 00 
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good correlation with the RNA-seq data, acute GR re- 
cruitment peaks previously identified by ChlP-seq in Dex- 
treated adipocytes [32] were enriched in Dex-induced 
clusters 1, 6, 7 and 9 (Figure 7 A, arrow up). Among Dex- 
regulated TFs, mouse genome-wide binding datasets are 
currently available for KLF4 (ChlP-seq), KLF2 (chip-on- 



chip), PPARy (ChlP-seq) and NFIL3 (chip-on-chip) [36-40]. 
Although KLF4 sites are enriched in the entire Dex/ 
LPS dataset compared to non-expressing genes, only in 
cluster 8 (Dex-repressed genes) the enrichment level at- 
tains significance (Figure 7A). Conversely, KLF4 binding 
sites are underrepresented in the LPS -induced cluster 2. 
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Figure 6 GR binds genomic sites associated with genes encoding Dex-induced TFs. (A) ChlP-seq data from BMMO (black; see also Additional 
file 3: Figure S4), adipocytes (crosshatched) and C2C12 cells (red) reveal multiple Dex-dependent GR binding peaks associated with Dex-induced genes 
encoding TFs. The numbers under the putative GR binding sites indicate primer pairs used in ChlP-qPCR experiments in BMMO (Additional file 2: 
Table S3). (B) ChlP-qPCR-assessed GR occupancy at putative GR binding sites indicated in (A) prior (U) and upon (D) a 40-min stimulation with 100 nM 
Dex. (C) Scanning ChIP of the Klf4 genomic region reveals a putative GR binding site. Primer positions and amplified products are shown as black 
rectangles over the gene schematics. Putative sites harboring sequences similar to GR binding sites are labeled as GREl, 2,3 with the binding sequence 
shown on top. The capital letters indicate nucleotides fitting a consensus GRE. A green rectangle overlapping the Klf4TSS indicates the position of a 
KLF4 binding site identified in [36]. Data for each site represent average of 2 or more independent experiments; error bars are standard deviation. 



Although the KLF2 chip-on-chip dataset available was 
relatively small, several KLF2 targets were regulated by 
Dex and contained GR binding sites, including Klf2 itself, 
Klf4, Klf9 and TgfbS [37]. 

Two mouse PPARy datasets, one from differentiated 
adipocytes and one from resting MO are currently 



available [38,39]. Consistent with the previously re- 
ported role of PPARy in repression of inflammatory 
genes, PPARy binding sites are overrepresented in 
several LPS -induced clusters (clusters 2, 3 and 5; 
Figure 7 A) in both datasets. In addition, in the MO 
dataset, PPARy binding sites were enriched in Dex- 
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Figure 7 The distribution of GR, KLF and PPARy binding sites at Dex/LPS-responsive genes in BMMO. (A) For each expression cluster 
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tine gene and 15 Kb of flanl<ing sequences on eacli side (blue) was compared to those in random group of non-expressors (RPKM < 1, red), 
low expressors (1 < RPKM < 10, green) and all expressing genes (RPKM > 1, orange) in BMMO. (B) A select group of Dex-responsive genes 
contains clusters of GR (black - macrophages (this study), crosshatched - adipocytes [32], red - C2C12 cells [33]), KLF4 (blue; [36]) and PPARy 
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induced and -repressed clusters, 7 and 8, respectively 
(Figure 7A). 

The only available genome-wide dataset of Nfil3 bind- 
ing was acquired in a neuronal cell line [40]. Among 
Nfil3-occupied genes identified in that study, only four 
genes overlap with Dex/LPS -regulated dataset reported 
here; however, one of them, Tsc22d3 (GILZ), is a well- 
characterized GR target. 

To identify genes that might be under combinatorial 
control by GR and another Dex-responsive TF, we 
searched for loci that contained GR, KLF and/or PPARy 
binding sites located close to each other within a gene +/- 
15 Kb. Intriguingly, several GR targets including Klf2, 
Klf9, Cited2 and Mt2 contained tight clusters of binding 



sites for GR, KLF(4) and PPARy (Figure 7B) suggesting 
that these TFs may interact functionally or physically at a 
subset of GR-regulated genes. 

LPS down regulates a unique class of genes encoding the 
C2H2-KRAB gene expression regulators 

Even a brief LPS treatment results in a marked downreg- 
ulation of a large number of genes confined to clusters 
10-12, GO overrepresentation analysis of LPS -repressed 
genes revealed that many of them participate in the 
regulation of nucleic acid metabolism, gene expression 
and transcription however, detailed information on spe- 
cific functions of many of these proteins is lacking. Inter- 
estingly, 33 proteins in these clusters contained tandem 
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Figure 8 Many LPS-downregulated genes encode Zn-finger proteins containing KRAB domains. (A) Domain overrepresentation analysis 


among LPS-suppressed genes {Clusters 10, 11 and 12) revealed a large number of Zn-finger proteins (COG 5048) that, in addition, contain KRAB, 


SCAN and BTB domains. The list of LPS-downregulated Zn-finger proteins is shown with KRAB, SCAN and BTB domains indicated as red, blue and 


green boxes, respectively. (B) The kinetics of LPS-mediated suppression of select genes encoding Zn-finger KRAB domain-containing proteins. 


BMMO were treated with LPS for indicated time and the expression of indicated genes was assessed by RT-qPCR. 





zinc-finger motifs (COG:5048, p = 3.49*10"^^ Figure 8A). 
We confirmed that the expression of 10 of these Zn-finger 
proteins is rapidly downregulated by LPS in MO, and re- 
mains suppressed for up to 6 h of treatment (Figure 8B). 
Further analysis of domain architecture revealed that in 
the majority of these proteins, tandem Zn-fingers co- 
occur with domains such as Kruppel-Associated Box 
(KRAB, Pfam01352, p = 1,029*10"^^), BTB (Pfam00651, 
p = 5.68*10"^) and SCAN (Pfam 02023, p = 2.1*10"^ 
Figure 8A). KRAB is a Tetrapoda-specific domain that 
defines one of the largest sub-families of Zn-finger pro- 
teins [41] which are involved in nucleic acid binding 
and regulation of gene expression. Although the spe- 
cific functions of the majority of KRAB proteins with 
respect to innate immunity are not well studied, in a 
few characterized cases KRAB proteins have been asso- 
ciated with transcriptional repression, establishing re- 
versible patterns of histone and DNA methylation and 
reversible heterochromatization [42-44]. 



Discussion 

Glucocorticoids- and LPS-regulated gene expression 
programs 

GR is a ubiquitous ligand-dependent TF capable of elicit- 
ing highly divergent transcription programs with up to a 
third of protein-coding genes differentially expressed fol- 
lowing a 24-h glucocorticoid treatment [7]. Establishing 
the hierarchy of regulatory events upon prolonged hormo- 
nal exposure in individual cell types is challenging, which 
complicates both accurate mechanistic predictions and 
clinical decisions. Multiple GR ligands have been designed 
in an attempt to create a highly specific compound that 
selectively regulates desired subsets of genes. Mechanistic 
analyses of these ligands usually focus on a specific group 
of disease- relevant genes and often involve long-term 
treatments, which obscure primary and transient re- 
sponses to GR activation by a plethora of secondary path- 
ways. In the context of inflammation, both immediate and 
delayed regulatory events are clinically relevant as they 
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reflect typical glucocorticoid treatment modalities. We 
reasoned that by analyzing early transcriptomes elicited 
by the inflammatory and glucocorticoid exposure in 
MO, a clinically relevant cell type, we will be able to 
isolate a key set of immediate GR targets responsible 
for the delayed gene expression patterns. Our results 
indicate that early glucocorticoid- and LPS -dependent 
changes establish a highly organized program of gene 
expression with distinct groups of genes following co- 
operative and antagonistic regulation. As expected from 
previous work [22,34,45] a large group of LPS-induced 
genes that included among others inflammatory cyto- 
kines, was rapidly downregulated by glucocorticoids. 
Another group encompassing glucocorticoid-induced 
genes, some of which encode TFs (Cited2, Nfil3, Jdp2) 
or signaling proteins (Duspl, Tsc22d3), are involved in 
curbing inflammatory signaling [46]. We identified several 
previously unreported glucocorticoid-induced genes whose 
products are involved in signaling (Ccll7, lUSra and 
Fzd4), regulation of transcription (Klf2, Jdp2, NcoaS 
and Tiparp) and mRNA stability (ZfandS). Several of 
these genes add to the arsenal of anti-inflammatory me- 
diators regulated by GR. For example, KLF2 interferes 
with API- and NFKB-mediated transcription of Tnf and 
several chemokines including Ccl2, 3 and 4 [47]. Fur- 
thermore, Klf2 haploinsufficiency in mice results in an 
exaggerated inflammatory response and more severe 
disease in arthritis models [48]. CCL17, another previ- 
ously unreported glucocorticoid-induced chemokine, is 
a marker and promoter of the polarization of alterna- 
tively activated' M2 MO, which are considered anti- 
inflammatory and mediate tissue repair and wound 
healing [49,50]. In addition to repressing cytokine gene 
transcription, glucocorticoids downregulate expression 
of several TFs including Atf3, Junb, Irfl, Bcl3, Tgifl, 
some {e.g., Rela, Nfl<b2, Myc, Etsl) in the context of LPS 
induction, and, unexpectedly, Pparg An enrichment in 
positive regulators of inflammation and cell proliferation 
among Dex-downregulated TFs is consistent with the 
anti-inflammatory and anti-proliferative effects of gluco- 
corticoids. The role of GR in repression of the Pparg 
gene in MO has not been previously reported, the effect 
might be indirect and mediated by a well-established 
GR target GILZ [51], which may also account for the 
delay (Figure 5A). Finally, we described a previously 
overlooked group of LPS -downregulated genes encod- 
ing proteins with the C2H2 Zinc-fingers adjacent to the 
KRAB domain. Despite being one of the largest TF 
family, KRAB proteins remain poorly characterized. 
Among those whose functions were described, several are 
involved in transcriptional regulation, RNA and DNA 
binding and splicing [42-44,52]. KRAB domains interact 
with a scaffolding co-repressor TRIM28 (KAPl, TIF1|3) 
which in turn binds the heterochromatin protein 1, 



chromatin remodeler MI2A and H3K9-specific methyl- 
transferase [53]. Indeed, some KRAB proteins reportedly 
repress transcription by heterochromatin spreading [52]. 
Interestingly, several KRAB proteins have been linked to 
NR actions [54,55]. The role of KRAB proteins in inflam- 
mation is essentially unknown; however, genomic studies 
indicate that inflammatory signaling increases accessibility 
of large sections of the genome [56]. It is tempting to 
speculate that a broad downregulation of proteins in- 
volved in heterochromatin maintenance and spreading 
serves to increase DNA accessibility and inflammatory 
gene transcription. 

The dynamic response to GR activation is consistent with 
feed forward logic 

Functional relationships between GR and its targets are 
often classified as "direct", that involve GR recruitment 
to genomic binding sites associated with regulated genes, 
and "indirect", whereby primary GR-regulated factors, 
rather than GR itself, are responsible for activation of 
the downstream targets. Thus, the activation of these 
secondary targets is often described as sequential or de- 
layed. Such a model, however, cannot explain many in- 
stances of non-monotonous expression dynamics (see 
Figure 5) and non-linear response to varying hormone 
concentration of many GRE-driven genes [57]. The large 
number of shared neighbors, overrepresentation of TFs 
and their high interconnectivity in GR regulatory net- 
works (Figure 2) are consistent with more intricate regu- 
latory modalities such as FFL. Variations in kinetic 
parameters for participating TFs, target gene structure 
and activation/repression thresholds often lead to para- 
doxical responses to stimulation of the master TF with 
profound functional implications. C-FFLs serve as delayed 
response organizers that detect the duration/strength of a 
signal that activates the initiating TF [11,12]. Interestingly, 
the dynamics of Fkbp5 induction by Dex, characterized by 
a substantial post-exposure delay followed by a robust ex- 
pression (Figure 5A), is reminiscent of the C-FFL in which 
the jointly regulated gene is activated by both the master 
and intermediate TFs [12]. Although additional experi- 
ments are required to establish the precise mode of Fkbp5 
regulation, this gene is a known direct GR target that re- 
cruits GR to several GREs (Additional file 3: Figure S5). 

Incoherent loops are responsible for negative and posi- 
tive pulse generation, accelerated response and fold 
change sensing [11,13,14]. Here, we observed that sev- 
eral GR target genes exhibit both positive (Klf2, Tiparp, 
Tgfl33 and Mt2) and negative (Tgifl, Junb and Bcl3) 
pulse-like dynamics consistent with the I-FFL. In keep- 
ing with the role of a potential master regulator, GR 
binds to the GREs in regulatory regions of many of these 
genes (Figure 6). Furthermore, using a system of ordin- 
ary differential equations which describe FFLs in the 
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"fold sensors" model [13], we showed that Klf2 expres- 
sion is consistent with that of a gene under joint control 
of GR and a strong GR- activated repressor (Figure 5C). 
Several GR-activated genes are either known transcrip- 
tion repressors {e.g,, Klf4, 9, Nfil3, Perl and Jdp2) or 
may downregulate gene expression by destabilizing RNA 
transcripts (ZfandS, [58]). Curiously, the expression dy- 
namics of Klf9 fits closely with the computational pre- 
diction of an intermediate repressor in the GR-R-Klf2 
I-FFL (Additional file 3: Figure S3). GR is recruited to 
the Klf9 and Klf2 GREs as early as 40 min of Dex treat- 
ment. Both Klf9 and Klf2 regulatory regions also contain 
functional GAGGCGTGG KLF sites ([36], Figure 7B) 
which can be occupied by various TFs of the KLF family 
[59] including KLF9. Finally, in KLF9-KO macrophages, 
the induction profile of Klf2 loses the early peak 
followed by a decrease and acquires monotonous kinet- 
ics (Figure 5D) strongly suggesting a collapse of the II- 
FFL to simple GR-dependent activation. Interestingly, 
KLF binding sites are overrepresented in glucocorticoid- 
regulated genes and are located near GREs in several bona 
fide GR target genes suggesting that these factors may co- 
regulate a number of GR targets. 

KLF proteins in inflammation 

Both KLF2 and KLF4 have been implicated in myeloid cell 
biology. KLF2 inhibits monocyte activation by inhibiting 
NFkB activity, which correlates with decreased expression 
of multiple cytokines and HIFla, a TF that regulates mye- 
loid cell response to bacterial infection and reactive oxygen 
species [60]. Consistent with the anti-inflammatory role of 
KLF2, mice hemizygous for Klf2 have elevated levels of in- 
flammatory mediators, such as CCL2 and PTGS2 (COX-2). 
By extension, in KLF2-deficient mice, the manifestations of 
both Me-BSA- and ILlp-experimentally induced inflamma- 
tory arthritis are more severe [48]. 

KLF4 is involved in inflammatory monocyte differenti- 
ation [61,62] and in MO polarization toward the M2 
anti-inflammatory phenotype [63]. KLF9 can act as 
either a transcriptional activator or a repressor [64,65], 
however its role, if any, in inflammation has not been 
described. We showed here that GR regulates Klf genes 
with distinct temporal dynamics and proposed that 
KLF9 may act as a GR-induced Klf2 repressor. Thus, it 
is tempting to speculate that GR anti-inflammatory ac- 
tivities rely in part on the activation of Klf genes whose 
products regulate transcription of additional targets in 
concert with GR. Indeed, glucocorticoids and KLF4 
regulate partially overlapping set of genes during epider- 
mal barrier establishment in embryogenesis [66]. The 
proximity of the GREs and KLF binding sites in the gen- 
ome suggests an intriguing possibility that GR and KLFs 
interact functionally or physically. Curiously, a func- 
tional interaction with the I-FFL logic has been reported 



for GR and another member of KLF family, KLF15 [67]. 
Although KLF15 is not expressed in MO, our studies 
strongly suggest extensive crosstalk between GR and 
other KLF family members in the innate immune cells. 

Conclusions 

Anti-inflammatory activities of glucocorticoids involve down- 
regulation of inflammatory mediators and activation of 
various anti-inflammatory genes. The early glucocorticoid- 
driven transcriptome in MO contains an unusually 
large number of genes coding for transcriptional regu- 
lators. Temporal dynamics of hormone-regulated gene 
expression is consistent with feed forward logic sug- 
gesting that GR and GR-induced TFs jointly regulate 
GR target genes. In particular, our data suggest that GR 
is rapidly recruited to and activates genes encoding sev- 
eral members of the KLF family of TFs with profound 
anti-inflammatory activities, such as Klf2 and Klf4. Fur- 
thermore, GR appears to regulate Klf2 expression via 
the GR-Klf9-Flf2 Il-FFL. We propose that by acting as 
a hub for highly branched regulatory networks and acti- 
vating genes encoding TFs to propagate the initial sig- 
nal, GR coordinates anti-inflammatory responses. 

Methods 

Mice and macrophage cultures 

C57BL/6 mice (NCI, Charles River Laboratories) were 
maintained in the Hospital for Special Surgery Animal 
Facility in fuU compliance with institutional guidelines 
approved by the HSS Animal Care and Use Committee. 
Klf9-KO mice [30] were generously provided by Dr. R. 
Simmen (U. of Arkansas). BMMO were prepared from 
8-12 weeks old mice as in [23]. For RNA-seq, BMMO 
from two independent mice were treated with vehicle, 
Dex (100 nM), LPS (10 ng/ml) or LPS + Dex for 1 h. For 
qPCR and ChIP analyses, BMMO were treated as above 
for time indicated in Figure Legends. 

RNA isolation, RT-qPCR and RNA-seq 

Total RNA was isolated using the RNeasy kit (Qiagen). 
0.25 (ig of total RNA was used for random primed 
cDNA synthesis which was performed with M-MuLV re- 
verse transcriptase (NEB) according to the manufac- 
turers recommendations. Quantitative PCR (qPCR) was 
performed using Maxima Sybr Green/ROX/ 2x master 
mix (Fermentas) on StepOne Plus real time PCR system 
(ABI) and analyzed using bbCt method as described pre- 
viously [22] with Hprt or Actl as a normalization control. 
Primer pairs are listed in Additional file 2: Table S3. RNA- 
seq is described in (Additional file 1). 

ChlP-qPCR and ChlP-seq 

BMMO were incubated -/+100 nM Dex for 40 min and 
ChlPs were performed as in [23] using N499 [22] and 
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scl004 (Santa Cruz Biotechnology) anti-GR antibodies, 
or normal rabbit IgG as a background control The data 
for each recruitment site was normalized to non-specific 
signals at the unrelated 28S ribosomal gene. Primer pairs 
are listed in (Additional file 2: Table S3). ChlP-seq is de- 
tailed in (Additional file 1). 

Gene association network construction, analysis and 
expression data modeling 

GeneMANIA algorithm was used to build gene associ- 
ation networks. Starting from a gene list of interest {e.g,, 
combined list of Dex up- and downregulated genes), 
geneMANIA algorithm creates a consensus network and 
predicts gene functions based on integration of multiple 
prebuilt gene association networks. The detailed descrip- 
tion of Gene association network construction, analysis 
and expression data modeling is in the Supplemental In- 
formation section. 

Additional files 
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